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I. INTRODUCTION 



The impact and spreading of liquid drops on solid substrates is the key element of a 
range of technological processes. Examples include spray cooling of surfaces, crop spraying, 



spray coating, solder jetting and DNA synthesis (Grissom and Wierum, 1981; Bergeron 



, 2000 




Maier et ai, 


2000) 



has been overwhelmingly devoted to the behaviour of millimetre-sized drops (Lesser and 
Field, 1983; Rein, |1993[ Yarin, 2006) where the spatio-temporal scales of interest allow 
experiments to be performed routinely (Rioboo et a/., 2002[ ). However, increasingly, there 
is an interest in the dynamics of microdrops, whose behaviour is critical to the functioning 



of a number of microfluidic devices (Squires and Quake, 2005). The ink-jet printer is one 



such device which, as well as being utilized in the graphic arts, is beginning to become 



a viable alternative to traditional fabrication methods ( |Gao and Sonin , 1994; Hong and 



Wagner, 1999; Calvert, 2001; Burns et a/., 2003), such as in the cost effective printing 



of P-OLED displays (Singh et a/., 2010) or the building of complex 3D structures through 
additive manufacturing (Derby, 2010). In such processes, the interaction of the microdrops 
with the solid substrate on which it impacts is directly related to the quality of the product, 
so that it is important to be able to predict and understand the behaviour of microdrops in 
such situations. 

Much research, currently confined to the millimetre scale, has focused on how the wet- 



tability of a solid can dramatically influence a drop's dynamics (Renardy et a/., 2003) and 
how surfaces can be designed with areas of high and low wettability in order to control the 
spreading of drops on them and enhance the precision of the process to correct for inherent 



inaccuracies in droplet deposition (Mock et a/., 2005). Given the even larger surface-to- 



volume ratio of microdrops, so that surface effects become more dominant, the possibilities 
for flow control using a pre-patterned solid substrate are significant; however, thus far very 
few experimental or theoretical papers have addressed this promising new area. 

The dynamics of microdrops as they impact on the solid substrate is difficult to observe 



experimentally (Dong et a/., 2007; van Dam and Clerc, 2004), especially with a sufficiently 



high temporal resolution, and characteristics such as the flow field inside the drop or stress 
acting on the substrate are almost impossible to measure. Therefore, it becomes necessary to 
rely on a theory which, once validated against the data from experiments on millimetre-sized 



drops and other relevant free-surface flows, would allow one to obtain reliable information 
about this process. 

So far, the main emphasis in research on microdrop spreading has been on the influence of 



additional physical effects on the drop's dynamics such as heat and mass transfer (Bhardwaj 



and Attinger, 2008); polymeric properties of the liquid (Perelaer et al, 2009); entrapment 



of bubbles under the drop (van Dam and Clerc, 2004); spreading on, and imbibition into, 



porous media (Clarke et a/., 2002; Holman et ah, 2002); evaporation (Lim et a/., 2009); 



and solidification (Attinger et al, 2000). The focus of our research programme is to begin 
by accurately capturing the process of dynamic wetting, which is the key physical effect 
in the drop impact phenomenon, and develop a benchmark numerical platform capable 
of incorporating complex mathematical models that describe the essential features of this 
process. Once this aspect is resolved, additional physical factors such as heat transfer, 
interaction with external fields, etc, can be built into the framework as required. An example 
of one such additional physical effect was given in |Sprittles and Shikhmurzaev! |2007[ [2009, 
where the model was extended to account for surfaces of inhomogeneous wettability. 

The issues surrounding the modelling of dynamic wetting flows are well known and have 



Dussan, 


1979 




Blake, 


2006 


) and monographs (e.g. 


de 



Gennes, 1985; Shikhmurzaev, 2007| ), with a general discussion of the merits of microscopic, 



mesoscopic and macroscopic modelling approaches given in Shi khmurzaev , 2011 and the 
discussion notes that followed. Here, we focus on the self-consistent framework of continuum 
mechanics where, in particular, it has been established that the classical model of fluid 



mechanics must be modified to allow for a solution to exist (Huh and Scriven, 1971; Dussan 



1979). A common way to achieve this is to use a 'slip model' where the no-slip condition 



is modified to allow for slip between the liquid and the solid near the moving contact line, 



e.g. using the Navier condition QNavier , 1823), whilst the contact angle is prescribed as a 



function of the contact-line speed and material constants (e.g. Greenspan, 1978). When 
incorporated into numerical software, such models have been shown to produce adequate 
results for the spreading of millimetre-sized drops at relatively low impact speeds, where 
experiments can be easily analyzed to allow for the development of a semi-empirical analysis 



, 1996 




Yokoi et ai, 


2009) 



An open question is whether the models that have been specifically developed for 
millimetre-sized drop dynamics can predict the behaviour of drops across a range of scales, 



i.e. towards micro/nanodrops, at a range of impact speeds. 

A step in the direction of answering this question was taken in the studies of |Bayer and 



Megaridis , 


2006 


and 


Sikalo et al , 


2005 



Their results show that, even with millimetre-sized 
drops, the contact angle is not simply a function of the contact-line speed, but is actually 
determined by the entire flow field. In other words, the dependence of the contact angle 
on the contact-line speed for given materials of the system is non-unique; it varies with the 
speed of impact, i.e. it depends on the particular flow. For an illustration of this effect we 



refer the reader to Figure 14 in Bayer and Megaridis, 2006 and Figure 9 in Sikalo et al. 



2005| This effect of the flow field on the contact angle is well known in the process of 
curtain coating where it has been termed the 'hydrodynamic assist of dynamic wetting' 



(Blake et a/., 1994; Clarke and Stattersfield, 2006). Similar dependencies of the dynamic 



contact angle on the flow field have been noted in the spreading of a liquid between parallel 



plates (Ngan and Dussan V, 1982), the imbibition of liquid into capillaries (Sobolev et al. 



2000, 2001b and in the coating of fibres (Simpkins and Kuck, 2003); but these flows are 



yet to undergo the level of scrutiny that the curtain coating process has where it has been 
shown using careful finite element simulations that the effect cannot be described using 



any interpretation of the conventional 'slip' models (Wilson et a/., 2006). Put simply, all 



currently available computational software, which implement the 'conventional' (i.e. slip) 
models, are unable to describe this key physical effect which has already been shown to be 
critical for the optimization of one industrial process and, as experiments suggest, is critical 
to the understanding of microdrop impact and spreading phenomena. 

Currently, the only model able to predict the influence of the flow field on the contact 
angle is the interface formation model. This model is based on the simple physical idea 
that dynamic wetting, as the very name suggests, is the process in which a fresh liquid-solid 
interface (a newly 'wetted' solid surface) is formed. The model is described in detail in 
Shikhmurzaev[ |2007| and has already been shown to be in agreement with experimental data 



, 1993, 


1997, 


1996, 


2005a|blc 



Sprittles and Shikhmurzaev, 2007, 2009). Notably, recent benchmark simulations for the 



technologically relevant phenomenon of a liquid-gas meniscus propagating through a cylin- 



drical capillary tube ( |Sprittles and Shikhmurzaev , 2012b) confirm that by modelling dynamic 
wetting as an interface formation process, we are able to demonstrate the effect of the flow 
field on the contact angle so that the aforementioned flows exhibiting assist-like behaviour 



can now be studied using our computational tool. 

A further advantage of the interface formation model over conventional approaches is 
that it is able to naturally account for the influence of variations in the wettability of the 
solid surface on the bulk flow. Specifically, even when there is no free surface present, 



the interface formation model predicts that a single (Sprittles and Shikhmurzaev, 2007) or 



multiple (Sprittles and Shikhmurzaev, 2009) changes in the wettability of the solid substrate 
can disrupt a shear flow parallel to that solid. The obtained results are in qualitative 
agreement with the predictions of molecular dynamics simulations (|Priezjev et a/., 2005 



Qian et a/., 2005) and will become more important as the scale of the flow is reduced, i.e. 



as we consider micro and nanofluidic flows. 

Computation of dynamic wetting flows is complex: besides the effects of capillarity, 
viscosity and inertia, one must also capture the physics of wetting which typically occurs on 
a length scale much smaller than that of the bulk flow. As explained in detail in |Sprittles and 



Shikhmurzaev[ 2012a[ the majority of publications in the field fail to accurately account for 



the wetting dynamics on the smaller length scale and, consequently, make it impossible to 
distinguish physical effects from numerical errors in their results. Such codes may provide 
realistically-looking results for millimetre-sized drops, where the accurate computation of 
the bulk dynamics may be sufficient, but on the micro/nano scale, where an increasing 
surface-to- volume ratio means that surface effects become more important, such codes are 
unreliable as the free surface shape near the contact line, specified by the contact angle, is 
not accurately determined. 

The first steps in the development of a benchmark code for such flows was described in 
great detail in |Sprittles and Shikhmurzaev] |2012a[ where a user-friendly step- by-step guide 
is provided to allow the reader to reproduce all results. This framework was extended in 
Sprittles and Shikhmurzaev , 2012b| to allow for the simulation of dynamic wetting as an 



interface formation process, with the approach robustly tested by checking its convergence 
both under mesh refinement, allowing practical recommendations on mesh design to be 
made, and, in limiting cases, to analytic asymptotic results. Furthermore, the developed 
computational tool was used to predict new physical effects, such as a new type of dependence 
of the dynamic contact angle on the flow geometry, and was seen to describe experimental 
data for the imbibition of a liquid into a capillary tube exceptionally well and significantly 



better than the Lucas- Washburn approach (Washburn, 1921). This situation is in complete 



contrast to commercial software which have not been validated as thoroughly for this class 



of flows and have been shown to converge to the wrong solution for similar flows (Hysing 



2007). 



In contrast to the flows considered in |Sprittles and Shikhmurzaev] |2012a|b[ impacting 
microdrops undergo extreme changes in shape. Therefore, having outlined the model in 
Sec. \U\ and the numerical approach in Sec. [TTTJ, in Sec. [lV]to verify the code's accuracy for 
large changes in free surface shape we compare our calculations to a known test-case from 
the literature of oscillating liquid drops, where reliable results exist for an unsteady problem 
in which inertia, capillarity and viscosity are all prominent; it is also a problem of significant 
interest in its own right. 

Having verified the code's ability to accurately simulate the dynamics of liquid drops, in 
Sec. [V] we demonstrate its capabilities for drop impact and spreading phenomena. Here, we 
focus on the main physical effects and, in particular, the influence of a substrate's wettability 
on an impacting drop's behaviour as well as the ability of our computational tool to extract 
information that is hidden to experimental analysis and which allows new paths of enquiry 
to be proposed. Most promisingly, in Sec.[VT]a novel mechanism of flow control is developed, 
based on pre-patterning a substrate with regions of high and low wettability, that allows the 
final shape of the drop to be controlled using only slight alterations in the impact speed. 



II. MODELLING OF DYNAMIC WETTING PHENOMENA 



Consider the flow of an incompressible Newtonian liquid, of constant density p and vis- 
cosity /x, surrounded by a dynamically passive gas of a constant pressure p g) so that the 
continuity and momentum balance equations are given by 

du 



Vu = 0, 



P 



dt 



+ u - Vu 



- Vp + /iVu + pg, 



i) 



where t is time, u and p are the liquid's velocity and pressure, and g is the gravitational 
force density. 

The interface formation equations, which act as boundary conditions for the bulk equa- 
tions ([I]), and have been described in great detail in Shikhmurzaev, 2007, are here simply 
listed with a very brief description given below. These equations consider interfaces as 
'surface phases' characterized by the surface density p s , surface velocity v 5 with which the 



6 



density is transported and the surface tension a which can be viewed as a 'two-dimensional 
pressure' taken with the opposite sign. On a liquid-gas free surface, we have 

df 



pn • 



dt 

Vu+(Vu) T • (I - nn) + Voi = 0, 

p(u-v s 1 )-n=(p s 1 -p s le )T-\ 



vf • V/ = 0, 



-p + pn- 



(2) 



Vu+(Vu) T -n = <7iV-n, (3) 



dt 



+ V-(p s y i ) = -(p{-p{ e )r- 1 , 



(4) 



4/3 (v*n - u,|) = (1 + 4af3) V^, a, = 7 (^ 0) - p>), (5) 

where the a- priori unknown free surface is f(r,t) = 0, with the inward normal n = > the 
metric tensor of the coordinate system is I and the subscript 1 1 denotes components parallel 
to the surface, which are obtained by convolution with the tensor (I — nn). Subscripts 1, 2 
refer to variables on the free surface and liquid-solid interface, respectively. At a stationary 
liquid-solid interface, the equations of interface formation have the form 



• n = 0, pn- 
p(u- v*) n = (p* -p s 2e )T~\ 

V 2|| ~ ^ U ll = «V<7 2 , 



Vu+(Vu) T • (I - nn) + |W5 = 



dp s 2 
dt 



+ v- (pM) = -0^-/4) r-\ 



(6) 



(7) 



02 = 7(P( ) - P\)- (8) 

Boundary conditions ([2])-(j8]) are themselves differential equations along the interfaces and 
are in need of boundary conditions at the contact line. At a contact line where a free surface 



meets a solid, we have continuity of surface mass flux and the Young equation (Young, 1805), 



which balances the tangential components of the forces due to surface tensions acting on 
the contact line and hence determines the dynamic contact angle 64: 



Pi ( v i|| - U c ) ■ ni! + p\ (v||| - U c ) • m 2 = 0, a 2 + a x cos 6 d = 0. 



(9) 



Here nij are the unit vectors normal to the contact line and inwardly tangential to surface 
i = 1,2, the velocity of the contact line is U c and the surface tension of the solid-gas interface 
is assumed to be negligible. 

On an axis of symmetry, with normal vector n a , the usual conditions of impermeability 
and zero tangential stress are applied 



u • n a = 0, 



n a • 



Vu+(Vu) T (l-n a n a ) = 0, 



(10) 



7 



together with conditions on the smoothness of the free surface where it meets such an axis 
n n a = and the absence of a surface mass source/sink for the interface formation equations 
vfj • n a = 0. 



On the free surface, the standard kinematic equation is ^ and the equations balancing 
the tangential and normal stress acting on the interface from the liquid and gas with the 
capillary pressure are On the liquid-solid interface, equations ^ state that the solid is 
impermeable and that the tangential component of the bulk velocity satisfies a generalized 
Navier condition which shows that slip on the liquid-facing side of the liquid-solid interface 
can be generated by both tangential stress on the interface and variations of surface tension 
in the interface. The model takes into account the mass exchange between the bulk and 
surface phases, in Q and ([7]), that are associated with the relaxation of an interface with 
surface density p s towards its equilibrium state p s = p s e on characteristic relaxation time r. 
The first equation in ^ and ^ shows that the tangential components vjj of the surface 
velocity are driven both by the bulk motion of the fluid and by gradients in surface tension 
along the interface. The second equation in ^ and ^ is the surface equation of state, 
which here is taken in its simplest linear form, that determines the surface tension along 
the interface from the surface density and hence allows one to find the contact angle from 
the Young equation in ([9]). Therefore, the mechanism by which assist occurs becomes clear, 
with the flow able to influence the value of the surface tensions at the contact line and hence 
alter the contact angle. Estimates for the phenomenological material constants a, /3, 7, r 
and have been obtained by comparing the theory to experiments in dynamic wetting, 
e.g. in |Blake and Shikhmurzaev] |2002 , 



Notably, the wettability of the solid substrate is controlled by the equilibrium surface 
density on it p s 2e) which is directly related to the surface tension in the liquid-solid inter- 
face via the equation of state ^ and hence to the equilibrium contact angle 9 e via the 
Young equation Then, as shown in Sprittles and Shikhmurzaev, 2007, variations in the 
wettability of the substrate, considered in Sec. [VTJ can easily be built into the model by 
prescribing a spatial dependence for p s 2e . 



III. A COMPUTATIONAL FRAMEWORK FOR DYNAMIC WETTING 
PHENOMENA 



A framework for the simulation of dynamic wetting flows as an interface formation pro- 
cess has been developed in |Sprittles and Shikhmurzaev] |2012b[ which builds upon previous 



work (Sprittles and Shikhmurzaev, 2012a) where the approach was formulated for the math- 



ematically less complex conventional models. These papers provide a step-by-step guide to 
the development of the code, curves for benchmark calculations and a demonstration of our 
platform's capabilities, by showing excellent agreement with experiments on capillary rise, 
so that, here, we shall just recapitulate the main details. 

The CFD code is based on the finite element method and uses an arbitrary Lagrangian- 
Eulerian mesh design (Kistler and Scriven, 1983; Heil, 2004; Wilson et al , 2006[ ) to allow the 
free surface to be accurately represented whilst bulk nodes remain free to move (snapshots 
of the mesh can be seen below in Figure [6]) . This mesh is based on the bipolar coordinate 
system and is graded so that exceptionally small elements are used near the contact line to 



accurately capture the physics of interface formation there (see Sprittles and Shikhmurzaev 



2012b for estimates on the smallest element size required for given parameters), whilst in 
the bulk of the liquid larger elements are adequate and ensure the resulting problem is 
computationally tractable. 

The result of our spatial discretization is a system of non-linear differential algebraic 



equations of index two (Lotstedt and Petzold, 1986) which are solved using the second- 



order backward differentiation formula, whose application to the Navier-Stokes equations is 
described in detail in |Gresho and Sani[ |1999[ using a time step which automatically adapts 
during a simulation to capture the appropriate temporal scale at that instant. For example, 
the time steps required in the initial stages after impact of the drop are small compared to 
those used in the later stages where the drop often gently oscillates around its equilibrium 
position on a much longer characteristic time scale. 

The described CFD code has been thoroughly validated for dynamic wetting phenomena 
in which the changes in free surface shape are moderate, but is here applied to microdrop 
impact and spreading where huge changes in free surface shape can be observed. Therefore, 
continuing in the spirit of our previous publications, the code is first validated for this class 
of flows which, additionally, allows us to demonstrate how easily the the framework can be 



9 



adapted to handle a variety of free surface flows. 



IV. OSCILLATING DROPS: VALIDATION OF THE UNSTEADY CODE 
FOR LARGE FREE-SURFACE DEFORMATION 

Consider the free oscillation of liquid drops in a parameter regime that will allow us 
to test our results against known numerical solutions for a high deformation unsteady flow 
which exhibits the competing forces of capillarity, inertia and viscosity. For small oscillations 



analytic results exist (Rayleigh, 1879), but for arbitrary viscosity and deformation numerical 



methods are required. The test case considered is for the standard model, which can be 
easily obtained from the interface formation model's equations for the free surface ([2])-([5]) 
by setting the ratio r/T of interfacial relaxation time r to bulk time scale T to zero (see 



Sprittles and Shikhmurzaev, 2012b for details), so that the same code can be used. In this 



case, the flow is fully characterized by the Reynolds, Stokes and capillary numbers, which 
are, respectively, 



Re 



pUL 



St 



pgL* 



Ca 



(11) 



II flU ' (Tie 

where L and U are the scales for length and velocity, and a\ e is the equilibrium surface 
tension on the free surface. Here, we consider oscillation in zero-gravity so that St = 0. 
Parameters are chosen to allow a comparison of our results with the numerical studies 



reported by other groups in (Basaran, 1992) and (Meradji et al, 2001). To do so, we run two 



simulations of the axisymmetric oscillation of a drop with Re = 10, 100 and Ca = 0.1,0.01 
(so that We = ReCa = 1 for both simulations), respectively. We consider the oscillation of 
a viscous liquid droplet whose initial shape is most naturally represented in spherical polar 
coordinates {R ) a ) Lp) ) with the origin located at the centre of the drop so that the drop 
surface is given by 

XG = /(a,t)e fl , (12) 

where e# is a unit vector in the radial direction. Then, as shown in Figure [T}i, a = /(0,i) 
is the length of the semi-major axis and b = /(7r/2, t) is the length of the semi-minor axis. 

In the benchmark test case, the drop is released from a shape whose deviation from 
a sphere is proportional to the 2nd spherical harmonic P 2 (cosa) = | (3 cos 2 a — 1), with 
coefficient of proportionality chosen to be 0.9, so that 



10 



Re=10 


— H - ^ 

Basaran 1992 
1 u 


Meradji et al. 




2001 


Present Work 


Ti 


2.660 


2.640 




2.656 


(o/6) Tl 


1.434 


1.432 




1.432 


Re=100 








Ti 


2.905 


2.930 




2.936 


(o/6)n 


2.331 


2.304 




2.305 



TABLE I. Comparison of current results with previous studies. 



/(a,0) =7[l + 0.9 P 2 (cosa)], (13) 

where 7 is a normalizing factor which ensures that the droplet has the correct non- 
dimensional volume, in our case 7 = An/ 3. 

We record the time 7\, and aspect ratio of the drop after one period (a/6)r 15 for Re = 
10, 100 and We = 1 in order to validate our code against the previous studies. For the 
results which are presented, a relatively crude mesh of 630 elements was used with a fixed 
time step of St = 0.001. Doubling the number of elements or reducing the time step by a 
factor of ten resulted in a change of less than 0.1% in both the time and amplitude recorded 



after one period. Significantly, the results in Basaran, 1992 were for 128 elements whilst 



Meradji et al] [2001] use an order of magnitude more elements. 

Figure [T] shows the evolution of the drop over one period for the Re = 100 case. The high 
deformation of the free surface is clear and one can observe that at the end of the first period 
the drop, whose equilibrium shape is a sphere, has its amplitude of oscillation damped by 
viscosity. Figure [2] shows the time-dependence of the aspect ratio a/6, which is recorded 
after one period t = T. It should be pointed out that the kinks in curve 1 of Figure [2] are 
not numerical artifacts; they are associated with the high deformation regime and can be 
seen in the previous studies. 

Our results in Table [I] are seen to be in good agreement with both studies. The values 



align most closely with those of |Meradji et al\ |2001[ which is reassuring given the greater 



mesh resolution associated with this study. 

Thus, it has been demonstrated that our numerical framework is able to provide accurate 
results for complex unsteady free-surface flows. The oscillation of liquid drops is a problem 

11 



(a)t=0 



(e)t=1.68 





(b)t=0.42 



(f)t=2.09 




(c)t=0.84 





(g)t=2.51 



6. 



(d)t=1.26 ( h)t = 2 - 93 
FIG. 1. Evolution of a liquid drop with Re=100 and We=l over one period (T = 2.93). 
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FIG. 2. Aspect ratio a/6 of two drops over a number of periods with curve 1 obtained using 
Re=100, We=l whilst curve 2 is for Re=10, We=l. 



of interest in its own right, and at this point we could look at comparing our results to 



experiments in the literature (Wang et a/., 1996), to probe newly proposed analytic models 



for decay rates (Smith, 2010) or to consider the influence that including interface formation 



physics may have when oscillations are of high frequency and the interface is forming and 
disappearing at a significant rate. All of these avenues of investigation are being pursued 
but lie beyond the scope of the present paper, and we now turn our attention to the code's 
capabilities at describing drop impact and spreading phenomena. 



V. MICRODROP IMPACT AND SPREADING 

To study the key physical effects in microdrop impact and spreading phenomena in a range 
applicable to inkjet printing technologies, we consider simulations of impact on hydrophilic, 
hydrophobic and patterned substrates. In particular, we see that our code is able to account 
for the two extreme outcomes of the drop impact and spreading process, i.e. deposition and 
rebound; to recover information about the drop's dynamics which is currently unobtainable 
from purely experimental analysis; and to predict new methods for flow control on chemically 
patterned surfaces. 

Consider a microdrop of water of radius L = 25 fim which impacts a solid substrate 
at U = 5 m s _1 , with the subsequent motion considered axisymmetric. Then, the non- 



13 



dimensional parameters are Re = 130, Ca = 0.07, St = 0.001 and estimates for the interface 



formation model's parameters are taken from |Blake and Shikhmurzaev[ |2002[ All that 
remains to be specified is the wettability of the solid substrate, which is characterized by 
the equilibrium contact angle 9 e that a free surface forms with the solid. 

Two simulations are shown in Figure [3] for the impact of a drop on substrates of different 
wettabilities, 9 e = 60° and 9 e = 130°, respectively, with an associated link to a movie 
here: |Movie[ The evolutions of the contact line radius r c and and apex height z a are 
given in Figure [4| Hereafter, for brevity we use the term 'apex' for the point located at 
the intersection of the axis of symmetry and the free surface; as Figure [7] shows, it isn't 
necessarily the highest point of the free surface. 

In Figure [3] it can be seen that in the early stages of spreading when inertia is dominant, 
roughly until t — 1 (with one dimensionless unit corresponding to 5 [is for the drop on which 
our non-dimensional parameters were based) , the shapes of the two drops are indistinguish- 
able. As can be seen in Figure [5| the contact line is forced outwards as fluid is pushed out 
radially from the centre of the drop which is being driven vertically into the substrate by 
inertial forces. This causes a toroidal rim of fluid to form near the contact line, with the 
pressure plot in Figure [5] clearly showing the formation of a disturbance, travelling from the 
contact line towards the apex, which separates the growing rim of fluid from the bulk of the 
drop. Eventually, the drop starts to feel the wettability of the solid on which it is spreading 
and, since in both cases inertia has carried the drop past its equilibrium position, the contact 
line starts to recede. Noticeably, as can be seen in both Figure [3] and Figure [4j although 
the wettability of the substrate eventually begins to alter the position of the contact line; 
this is not felt by the apex until a much later time. In fact, the initial fall of the apexes are 
very similar, and it is only upon their recoil, around t — 4, that their paths begins to differ. 
When the drops begin to recoil, their motions differ quite significantly; notably, for the drop 
on the hydrophobic substrate, where the dewetting of the substrate occurs so quickly that 
the drop rebounds back off the substrate. By comparing the images at t — 3 and t = 4 
in Figure [3| one can see that the rebound is preceded by a jet emanating from the apex 
region which pulls fluid radially inwards towards the axis of symmetry. This second stage 
of spreading is seen to be on a much longer time scale than the initial stages after impact. 
The simulation has to be terminated as the drop is about to leave the substrate; extending 
the numerical platform to account for such behaviour is certainly viable. It is of interest 
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FIG. 3. Microdrop impact and spreading simulations at Re = 130, Ca = 0.07, St = 0.001. The grey 
substrate (left) is hydrophilic (0 e — 60°) whilst the green one (right) is hydrophobic (9 e — 130°). 
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FIG. 4. Evolution of the drop's contact line position r c (1, 2) and apex height z a (la, 2a) as a 
function of time for simulations at Re = 130, Ca = 0.07, St = 0.001. Curves 1 & la: 9 e = 60° and 
curves 2 & 2a: 9 e = 130°. 

to see that the drop's final shape is pear shaped and, indeed, this shape has been observed 



experimentally (Mao et a/., 1997). 

Snapshots of the mesh used during the computation of the obtained results are given 
in Figure [HJ where a relatively crude mesh is shown, allowing the elements to be easily 
identified. One can see that the mesh remains regular throughout extreme changes in free 
surface shape, like at t = 3, when the apex is very close to touching the substrate, and at 
t = 10.5, when the drop is close to leaving the substrate. 

Some of our computational framework's advantages over a purely experimental analysis 
of the phenomenon are as follows. First, it can recover information which is inaccessible to 
experiments; second, one can efficiently map the influence of the system's parameters on the 
drop's dynamics, and, third, it is easy to attempt new things without the cost of full scale 
laboratory experiments. 

As an illustration of the first of these advantages, as shown in Figure [5] and Figure [7| in 
our simulations we are able to see the entire shape of the microdrop for the whole simulation, 
whereas experimental images on microdrops are unable to show the dynamics of the apex 
as it disappears below the rim of fluid which surrounds it, as well as features experimentally 
unobtainable at these scales such as the flow field and pressure distribution inside the drop. 
It can be seen that the apex gets extremely close to touching the solid substrate, i.e. to 
dewetting the centre of the drop. This can also be seen from curve 2a in Figure [4j The apex 
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FIG. 5. Microdrop impact and spreading simulation at Re = 130, Ca = 0.07, St = 0.001 on a 
hydrophilic substrate (6 e = 60°). Left: Velocity vectors, Right: Pressure field. 



manages to recover just in time: as the contact line is receding, the apex re-emerges out of 
the centre of the drop in a jet-like protrusion (see Figure [3]). One could envisage that if the 
drop's apex does dewet the centre of the drop, then this additional dissipation of energy may 
inhibit the rebound of the drop off the hydrophobic surface, although, initial indications 
in the relatively narrow parameter space investigated thus far show that microdrops are 
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0.5 1 1.5 1 2 3 1 2 3 

FIG. 6. Snapshots of a relatively coarse computational mesh during the impact, spreading and 
rebound of the drop in Figure [3] off a hydrophobic substrate at times t = 3,6,10.5. Enhanced 
online. 



FIG. 7. Shape of the microdrops at t = 3 highlighting dynamics which are inaccessible to experi- 
ments. 

remarkably resilient to this dewetting feature. Understanding, and hence being able to 
control, this feature would be of significant benefit to processes which are constrained to 
using hydrophobic surfaces but would like to inhibit rebound. This aspect of the drop impact 
phenomenon will be the subject of future investigation. 

With regard to the second advantage of reliable numerical simulations of microdrops over 
experiment, determining for what parameter values a drop will rebound is an important 
piece of information, particularly when the substrate to be used in a given process has to 
be hydrophobic, and it is difficult to find this out from experiments where one cannot vary 
material parameters of the system independently. With our numerical tool, parameter space 
can be mapped quickly and efficiently, so that one can ensure the drop deposits on a given 
substrate by, say, artificially changing the viscosity of the liquid or reducing the impact 
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speed of the drop. Changing the substrate for the same microdrop impact parameters, we 
see that, for example, at 9 e = 110° no rebound of the drop is observed. 

Next, as an example of the code's cost-effectiveness with regard to the process, we consider 
the impact of drops on a chemically patterned surface which is custom built to enhance pre- 
determinable flow control on the drop's dynamics. 



VI. MICRODROP IMPACT AND SPREADING ON CUSTOM BUILT 
CHEMICALLY-PATTERNED SURFACES 



Consider how topologically different patterns of wettability on a substrate can be used 
to gain a required level of flow control on a drop once it has been deposited. Here, one 
such pattern is considered with the aim of changing the final shape of a drop for the same 
liquid-solid combination by only slightly changing the impact speed, an outcome which is 
impossible on a homogeneous substrate. To do so, we pattern an otherwise wettable solid 
[6 e = 60°) with a circle of nonwettable substrate (0 e = 110°) of radius 1.52 times that of the 
drop's initial radius (Figure [8]). 

From Figure [8] and Figure |9j and the movie of the simulation whose link is here Movie 



we see how the desired flow control becomes realizable. On the patterned surface, the 
equilibrium radius of the area wetted by a drop impacting at 4 m s _1 is found to be r c = 1.03, 
whilst for a 5 m s _1 impact it is r c = 1.61. This occurs because for the lower impacting speed 
the drop is unable to reach the edge of the hydrophobic area to access the more wettable 
region, and hence behaves as if it is on a homogeneous surface with wettability defined 
by 9 e = 110°. For the higher speed of impact the drop is able to reach the edge of the 
hydrophobic disc, and encounter the more wettable substrate, as can be seen by looking at 
curve 2 in Figure [9] at t = 2. This results in an increase in the wetting speed and causes the 
contact line to advance further. This is no guarantee that the drop's contact line will remain 
on the more wettable surface as the contact line could return to the hydrophobic solid, which 
in turn would enhance the dewetting process. However, from curve 2 in Figure |9| we can see 
that the contact line's recoil is relatively shallow, and, in the case considered, it approaches 
an equilibrium position without ever encountering the hydrophobic disc again. 

Thus, it has been shown that custom built substrates can be designed very quickly using 
our code to determine the specific details, such as the radius of the nonwettable inner circle, 
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(c)t=2 (f)t=15 

FIG. 8. Evolution of two drops impacting a patterned surface at (left) U = 4 m s _1 and (right) 
U = 5 m s _1 . The hydrophobic surface patterning (r < 1.52, 9 e — 110°) is green whilst the 
hydrophilic one is grey (r > 1.52, 9 e — 60°). 

that are critical to the success of the product. Such a computational tool could be used as a 
precursor to full-scale laboratory experiments and would help to vastly narrow the bounds 
on potential parameter regimes for a given requirement. Although we were able to choose 
different final contact-line radii by varying the impact speed, once the solid had been chosen, 
we had no control over what this radii would be. With an annulus of hydrophobic surface 
one would have the possibility of using the impact speed to control the final wetted area. 
However, the main new physical effect is the one considered here. 

VII. CONCLUSION 

The ability of the computational framework developed in |Sprittles and Shikhmurzaev , 
2012a, b, which for the first time models dynamic wetting as an interface formation process, 
to provide high-accuracy benchmark simulations for flows with large changes in free surface 
shape has been demonstrated. An initial study into the key physical effects of microdrop 
impact and spreading, applicable to the inkjet printing regime, has been performed, and the 
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FIG. 9. Radius of the contact line (1, 2) and apex height (la, 2a) as a function of time for two 
drops impacting a patterned surface at 1: U = 4 m s _1 and 2\ U = 5 m s _1 , respectively. The 
boundary between the surface characterized by 6 e = 110° (r < 1.52) and that defined by 6 e = 60° 
(r > 1.52) is marked with a dashed line. 

effect of the wettability of a substrate shown to be critical to the drop's dynamics, which 
can now be recovered in regimes hidden to experimental analysis. We have shown that the 
strong influence of the substrate's wettability on the drop's dynamics can be utilized to 
design chemically patterned surfaces allowing one to change considerably the final shape of 
the drop by only slightly altering the impact speed. This is an entirely new physical effect 
which deserves further investigation and fine tuning. 

In this article, our focus has been on showing the capabilities of our computational 
platform and highlighting some of the key physical effects of microdrop impact and spreading. 
Future work will be concerned with comparing our results to both existing experimental data 
in the literature, with the initial simulations showing excellent agreement with those observed 
Dong et~aT\ |2007[ and, in parallel, using our code to design theory-driven experiments. 



m 



The latter could identify parameter regimes in which bifurcations in the flow behaviour 
occur, such as rebound of the drop, dewetting of the centre, etc, and highlight where the 
differences between the predictions of models proposed in the literature for dynamic wetting 
will be most prominent. 

The finite element framework developed has already been shown to possess reasonable 



flexibility: it was used to consider flow over patterned surfaces (Sprittles and Shikhmurzaev 



2007, 2009), to simulate two-phase flow through a capillary (Sprittles and Shikhmurzaev 
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2012a,b) and here it was used to study the dynamics of free drops and their interaction with 



solid surfaces of varying types. This flexibility allows our research programme to branch 
out and simulate a whole array of different dynamic wetting flows, e.g. the coating of fibres 



(Ravinutala and Polymeropoulos, 2002), where the high speeds of coating in confined areas 



suggest 'hydrodynamic resist' (Sprittles and Shikhmurzaev, 2012b) may be present, or it 
can be applied to entirely different flows where interface formation has also been shown to 
be critical, such as the coalescence of liquid bodies, breakup of liquid jets, disintegration of 



thin films, and other phenomena (Shikhmurzaev, 2007). 
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